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Positive and negative Lyapunov exponents for a dilute, random, 
■ two-dimensional Lorentz gas in an applied field, E, in a steady 

^ ■ state at constant energy are computed to order E 2 . The results 

are: X± = A± — a±(qE/mv) 2 tQ where are the exponents for 
the field-free Lorentz gas, a + = 11/48, a_ = 7/48, to is the mean 
' free time between collisions, q is the charge, m the mass and v is 

the speed of the particle. The calculation is based on an extended 
Boltzmann equation in which a radius of curvature, characterizing 
the separation of two nearby trajectories, is one of the variables 
in the distribution function. The analytical results are in excellent 
' agreement with computer simulations. These simulations provide 

'-^ ■ additional evidence for logarithmic terms in the density expansion 

of the diffusion coefficient. 

^ ■ 

One of the outstanding problems in transport theory is to relate 
macroscopic quantities such as transport coefficients to microscopic quantities 
that characterize the chaotic dynamics of the system. This chaotic dynamics is 
responsible for the stochastic-like behavior which leads to normal hydrodynamic 
processes taking place in fluids. One approach is to consider the system as 
a Hamiltonian system, obeying classical mechanics, and to apply the escape- 
rate formalism which can be used to express transport coefficients in terms of 
Lyapunov exponents and Kolmogorov-Sinai entropies for trajectories in phase 
space on an appropriate fractal repelier [p. Another approach to this problem is 
to consider driven, thermostatted systems, in a nonequilibrium steady state where 
the system is subjected to an applied force as well as to a "Gaussian thermostat", 
which allows the system to reach a steady state by removing the heat produced 
by the applied force 0|§. Here one can relate the transport coefficients to the 
change in the sum of all, or in some cases, a pair , of the Lyapunov exponents 
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in the steady state. 

The system considered here is a dilute, thermostatted two-dimensional 
Lorentz gas, where a particle with charge q and mass m moves in an infinite 
random array of fixed hard disk scatterers of radius a and density n,na 2 <C 1, 
subject to a constant, uniform, external field E as well as to a frictional thermostat 
which maintains the speed of the particle at a constant value v. The system 
approaches a non-equilibrium steady state, characterized by an attractor in phase 
space. In this steady state, the following relation holds for the macroscopic 
diffusion coefficient, D, for the particle (or, equivalently, the conductivity a 
related to D by an Einstein formula D = m(v/q) 2 a ) and the Lyapunov exponents 
A+(e) and A_(e) §,§: D = - \im £ ^ [v 2 (\ + (e) + A_(e))]£" 2 , where e = (qE/vm). 
Here we report the first theoretical calculation of both A+(e) and A_(e), for 
this system, using kinetic theory. We also present a comparison with extensive 
computer simulations, as well as with the general relation, given above. We apply 
a method based on the extended Lorentz-Boltzmann (LB) equation developed by 
van Beijeren and Dorfman for computing the Lyapunov exponents [|5]||. 

The motion of the particle is described by a position vector r, and velocity v. 
Between collisions with the disks, the (non-Hamiltonian) equations of motion of 
the particle in the field and with an energy conserving Gaussian thermostat are 

x = v x = v cos 6; y = v y = v sin 9 

v x = ve- av x ; v y = -av y . (1) 

Here 9 is the angle the particle's velocity makes with the applied field, in the in- 
direction, and a = v x e/v represents the strength of the frictional force provided 
by the Gaussian thermostat 0, determined by the condition that the kinetic 
energy remain constant. 

The instantaneous change in velocity of the moving particle upon collisions is 
v' — v = — 2k(v ■ k), where v' is its velocity after collision, and k is the unit vector 
in the direction from the center of the scatterer to the point of impact. 

To determine the positive Lyapunov exponent, A + we first consider the 
separation of two diverging trajectories that start simultaneously at the same 
initial spatial point, but have slightly different initial velocity directions, specified 
by angles 9 and 9', with 9 (t = 0) = 9 , and & (t = 0) = 9 + 69 . We choose to 
measure the separation of trajectories along a line which is at all times t > 
perpendicular to the reference trajectory. Because the trajectories are curved by 
the thermostatted field, the intersection of this line with the adjacent trajectory 
does not occur at the point that a particle following the adjacent trajectory will 
reach at time t, even to first order in the infinitesimal quantities, such as 59 . 
In order to determine the Lyapunov exponents, one must take effects produced 
by the curving of the trajectories into account (see Fig. 1). The separation 
between the two trajectories at time t, SS (t), can be written as the product of 
a radius of curvature p(t) and S'9 (t), the difference between the velocity angles 
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of the two trajectories at the points of intersection with the perpendicular line, 
or pit) = S(t)/5'9 (t) . The radius of curvature changes continuously between 
collisions, and instantaneously at the collisions of the particle with the scatterers. 
Simple geometric considerations yield an equation for the rate of change of pit) 
between collisions as 

2 2 

p = v + (pe) cos 6 + sin 2 6. (2) 
The instantaneous change in the radius of curvature due to a collision is given by 

m 

112 

— = — + t, (3) 

p + p_ a cos (p 

where v ■ k = — cos 0, <fi is the angle of incidence at collision, and p_ and p + are 
the radii of curvature before and after collision, respectively. 

We now compute the positive Lyapunov exponent A + from the rate of 
separation of diverging trajectories. The result of Sinai [p]] for the positive 
Lyapunov exponent in a field-free Lorentz gas can be straightforwardly 
generalized to include the case under discussion: 

1 fto+T v I v \ 

\ + = Hm 1 / -y- dt= (-) . (4) 

+ T-^TJta P (t) \p/ ss {) 

Here we assume ergodicity and calculate the ensemble average in the steady state, 
denoted by the angular brackets with subscripts s.s., using a non-equilibrium 
steady state distribution function for the moving particle in the constant 
thermostatted applied field. The main idea is to consider an ensemble of similarly 
prepared systems, assume that the distribution function for the moving particle, 
in this ensemble, reaches a spatially homogeneous steady state (since there is 
no way to distinguish one spatial point from another in the ensemble average), 
and to derive and solve a Lorentz-Boltzmann (LB) equation for the distribution 
function for the moving particle f(v,p), where the variables include both the 
velocity and the radius of curvature describing the separation of trajectories of 
the moving particle and an adjacent trajectory, as described above. 
Thus we write 

\ + = dv dp-f(v,p) , (5) 
J Jo p 

assuming that / has been normalized to unity. The LB equation for / is 
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Since the dynamics between collisions is not Hamiltonian, the usual form of the 
streaming terms on the left hand side of the LB equation must be replaced by 
the form that reflects the total conservation of particles under the actual, non- 
Hamiltonian dynamics. The left hand side of the LB equation can be obtained 
from the equations of motion, Eq. ([TJ), and Eq. (g) for p 



v + pe cos 9 H j / 




(7) 

coll 



where we used the constant speed of the particle to denote f(v,p) by f(9,p). 
The right hand side of the LB equation is the change in / due to collisions given 
previously || as 

r' 2 aa, a,Ta<A (acos0)/2 \ r 

nav d<pcos<p dp 5 p - — — f{y,p) 

J-n/2 Jo \ 1 + (acos0)/2p'y 

~2navf(v,p). (8) 




coll 



To leading order in the density, the collision terms can be simplified further by 
approximating the factor 1 + (acos</>)/2// inside the 5 function by unity since p' 
is typically on the order of the mean free path, so that p' ^> a. Then Eq. (|7|) can 
be solved by expanding / as a power series in e and inserting the solution into 
Eq. (|). We obtain 

\ + (e) = \ -^t s 2 + O(e 4 ), (9) 

where to — ^l v i with i = (2na)~ l the mean free path of the moving particle, and 
Ao is given in Ref. ||. 

The calculation of A_(e) requires the study of trajectories that asymptotically 
converge to the reference trajectory of the moving particle rather than diverge 
from it. That is, the negative Lyapunov exponent can only be determined if one 
can find trajectories that lie on the stable manifold of the reference trajectory, 
which is typically difficult since almost all of the adjacent trajectories will 
eventually diverge from it. To overcome this difficulty we consider the time 
reversed motion of the moving particle ||. This allows us to consider diverging 
trajectories again. We thus consider the steady state distribution function f-(v, p) 
for trajectories of the particle with velocity v and radius of curvature p on the 
time-reversal of the stable manifold. The equation for /_ takes the form of an 
anti-Lorentz-Boltzmann (ALB) equation. This unusual form is dictated by the 
observation that if the moving particle and the scatterer with which it collides 
are uncorrelated before collision in the forward motion, then in the time reversed 
motion, the moving particle and scatterer will be uncorrelated after the collision. 
That is, to obtain the ALB, one must use the Stosszahlansatz for the exiting 
collision cylinders, rather than for those before the collision |||9] . The ALB then 
reads 
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v + pe cos 9 



p 2 e 2 



sin 2 6 /_ 



'df- 
dt 



(10) 
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where 



'df- 
dt 
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/ dp' 5(p - (a cos 0) 1 2) /_ (tf , p) 



J °°dp'f-(v,p>) 



The first term on the r.h.s. of Eq. (11) is the gain term. It is constructed by: (a) 
using the Stosszahlansatz for the collision cylinders centered on particles which are 
produced in collisions with velocity v, (b) requiring that the radius of curvature 
after collision be p, and (c) noting that before collision, the radius of curvature 
will typically be on the order of a mean free path, so that after collision p will 
be very close to (a/2) cos 8. The loss term in Eq. ( |TT1) can be obtained by noting 
that: (a) The rate at which particles with velocity v disappear due to collisions 
with angle of incidence is nav cos0 J °° dp'f_(v ', p'), and (b) the fraction of 
those which disappear having radius of curvature p is f-(v, p)[/o°° dp'f-(v, p')]~ l - 
These considerations lead directly to the r.h.s. of Eq. fllTD, and a more detailed 
explanation will be given elsewhere ||. 

Eq. (|10"D can be solved by imposing the requirement that 



dpf- (v,p) = 



(12) 



with f-(v) the time reversed steady state solution of the spatially homogeneous 
Lorentz-Boltzmann equation. To understand this condition, we note that when 
Eq. ( [TOD is integrated over all values of p to obtain an equation for the velocity 
distribution function /_ (v) one obtains an ALB equation with a collision operator 
that has the opposite sign from the usual Lorentz-Boltzmann equation, due 
to the fact that we are considering the time-reversed motion, and using the 
Stosszahlansatz after collisions, 0. For consistency, we then require that the 
steady state solutions of the p integrated LB and ALB equations be related by a 
simple time-reversal operation, as indicated by Eq. fll~2|) . 
The negative Lyapunov exponent is obtained from /_ by 



dv I dp-]- (v, p) 
o p 



(13) 



with /_ normalized to unity. The negative sign on the r.h.s. of Eq. ( pT3| ) is the 
result of the time reversal procedure. Again, the equation for /_ can be solved 
by expanding /_ in powers of e. We finally obtain 



-Ar 



(14) 
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Using the relation between the diffusion coefficient D and the sum of the 
Lyapunov exponents discussed earlier, we recover the correct low density value 
D = (3/8)t v 2 . 

It is important to compare these results with those obtained from computer 
simulations of the same system. For this purpose we distributed 10 5 non- 
overlapping scatterers randomly in a square simulation cell with periodic 
boundaries. Between collisions the trajectory of the moving particle is computed 



from an analytical solution |10[ of the thermostatted equations of motion. 
The collision points with the scatterers are determined numerically with an 
accuracy of 1CT 12 . The Lyapunov exponents were computed in tangent space 
by a simultaneous integration of the linearized equations of motion for the 
intercollisional streaming and an exact linearization of the map which relates 



the separation of trajectories after collision to that before collision JTTJ. In Fig. 



2 we show the deviation of the non-vanishing Lyapunov exponents from their 
equilibrium values as a function of the squared applied field. All numbers are made 
dimensionless by using a, v, m and q as the respective units for length velocity, 
mass, and charge. Two reduced scatterer densities n* = na 2 were considered, 
0.001 (circles) and 0.002 (squares). Each point is obtained by averaging over 
10 simulation runs with different scatterer configurations and a total of 4 x 10 6 
collisions for each run. The standard deviation for the exponents is ~ 0.1%. The 
lines refer to the theoretical expressions, Eq. (^) and Eq. (II). Both A + and A_ 
exhibit the predicted quadratic weak field behavior. These results are clearly 
consistent with the theoretical predictions for the field-dependent Lyapunov 
exponents, though at densities somewhat higher than those for which the LB 
results hold without density corrections. Fig. 3 shows D/(£v) as a function of 
e 2 £a/v 2 for n* = 0.001 (full circles) and 0.002 (open circles). The diffusion 
coefficient D is obtained from the conductivity o through the Einstein formula, 
and a is numerically computed from a = (qv x ) s , s ,/E, the ratio of the time- 
averaged dissipative current to the applied field. Alternatively, D can be obtained 
also from the Lyapunov exponents as indicated above. On the scale of Fig. 3 both 
methods yield indistinguishable results. 

We mention here that to obtain agreement between the simulation data for D 
and theory it was necessary to add to the low-density Lorentz-Boltzmann result 
(3/8)£v, the dotted top horizontal line in Fig. 3, known [[CJ moderate density 
corrections, which include essential contributions from logarithmic terms of the 
form n* In n*+0(n*). The numerical nonequilibrium results converge well to these 
corrected values for vanishing field indicated by the dashed in* = 0.001) and solid 
(n* = 0.002) horizontal lines. Thus the computer results confirm the presence of 
these non-analytic terms in the density expansion. 

As a check we computed the Lyapunov exponents for still smaller densities, 
n* = 0.0001, and fields < 10~ 5 with a very efficient direct-simulation Monte- 
Carlo method and found perfect agreement with our other simulation results 
quoted above. 
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We conclude by remarking that this work illustrates the power of kinetic 
theory methods for computing quantities of interest to both statistical mechanics 
and dynamical systems theory [p|,|6|,|I3|| . 
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Figure captions 



Fig. 1: Geometry for the curvature p{t) in the applied field. 

Fig. 2: Field dependence of the deviation of the Lyapunov exponents from 
their equilibrium value for the reduced densities n* = 0.001 and n* = 0.002. 
e = qE/vm, where E is the applied field, i is the mean free path, and q, m, and v 
are the charge, mass, and velocity of the particle, respectively, a is the diameter 
of the scatterers. 

Fig. 3: Dependence of the diffusion coefficient D on the field for the two 
densities n* = 0.001 and 0.002. All other quantities are as in Fig. 2. The horizontal 
lines are explained in the text. 
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